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Abstract 


This  report  presents  a  suboptlmal  boundary  estimation  algorithm 
for  noisy  images  which  is  based  upon  an  optimal  maximum  likelihood  problem 
formulation.  Both  the  maximum  likelihood  formulation  and  the  resulting 
algorithm  are  described  in  detail,  and  computational  results  are  given. 

In  addition,  the  potential  power  of  the  likelihood  formulation  is  demon¬ 
strated  through  the  presentation  of  three  simple  but  insightful  analyses 
of  algorithm  performance. 


Keywords 

Image  boimdary  estimation,  likelihood  maximization,  Markov  processes, 
tree  searching  algorithms. 


IMPLEMENTATION,  INTERPRETATION  AND  ANALYSIS  OF  A 
SUBOPTIMAL  BOUNDARY  FINDING  ALGORITHM 
1.  Introduction 

In  this  paper  we  present  an  algorithm  for  estimating  object  boundaries 
in  noisy  black  and  white  images.  The  algorithm  is  based  upon  a  maximum 
likelihood  Markov  process  state  estimation  formulation  developed  by  Cooper 
[1],  [2].  The  Images  treated  consist  of  a  constant  grey  level  object 
surrounded  by  a  constant  grey  level  background.  The  entire  picture  is 
assumed  corrupted  by  an  additive  white  Gaussian  noise  field.  The  algorithm, 
which  sequentially  maximizes  a  suboptimal  likelihood  function  to  obtain  a 
boundary  estimate  consisting  of  a  sequence  of  pixel  edge  elements,  serves 
to  Illustrate  the  tradeoffs  between  computational  feasibility  and  theoreti¬ 
cal  optimality.  As  la  demonstrated,  the  likelihood  formulation  also  pro¬ 
vides  a  framework  for  analysis  of  algorithm  performance  and  for  comparison 
with  other  algorithms. 

Our  algorithm  was  originally  motivated  by  a  similar  heuristic  tree 
searching  algorithm  introduced  by  Martelli  [3].  This  algorithm  estimates 
boundaries  by  minimizing  a  cost  function  containing  gradient  and  curvature 
components.  It  will  be  shown  that  this  algorithm  which  is  not  based  upon 
any  optimal  problem  formulation  has  potentially  poorer  performance  charac¬ 
teristics.  Actually,  our  approach  is  more  analogous  to  estimation  pro¬ 
cedures  used  in  control,  communication  and  information  theoretic  applications 
In  particular,  we  develop  a  Markovian  boundary  generation  model,  and  view 
boundary  estimation  as  a  problem  of  estimating  the  state  of  this  model  from 
noise  corrupted  measurements. 

As  discussed  by  Forney  [4],  classical  problems  such  as  convolutional 
coding,  frequency  shift  keying,  and  text  recognition  can  be  formulated 


In  the  Markov  state  estimation  framework.  Scharf  and  Elliott  [5]  have  also 
discussed  a  number  of  problems  In  signal  and  Image  processing  which  can 
be  solved  using  such  techniques.  A  conmon  feature  of  these  algorithms  Is 
the  formulation  of  a  likelihood  function  which  can  be  recursively  defined 
and  hence  sequentially  maximized.  Furthermore  maximization  of  this  likeli¬ 
hood  leads  to  a  maximum  a  posterioiH,  probability  (MAP)  estimate  of  the 
state  sequence.  In  many  cases  elegant  dynamic  programming  (sometimes  re¬ 
ferred  to  as  the  Vlterbl  algorithm  In  the  Information  theory  literature) 
algorithms  exist  for  the  sequential  maximization.  In  other  cases  where  the 
state  space  Is  large  other  graph  searching  algorithms  are  necessary.  As 
we  show  below,  although  It  is  a  simple  matter  to  derive  a  likelihood  func¬ 
tion  for  MAP  boundary  estimation.  It  Is  not  possible  to  calculate  this 
likelihood  recursively.  The  reason  for  this  Is  the  MAP  formulation  re¬ 
quires  use  of  all  the  picture  data.  When  object  boundaries  are  estimated 
sequentially  one  cannot  classify  all  picture  based  upon  a  partial  boundary. 
As  a  consequence  we  derive  a  suboptlmal  recursive  likelihood  function.  It 
Is  based  upon  the  optimal,  but  It  makes  use  of  only  picture  data  in  a  swath 
about  the  boundary.  Furthermore  due  to  the  size  of  the  state  space  we  are 
forced  to  use  a  suboptimum  modified  A*  (or  branch  and  bound)  graph  search 
algorithm  to  perform  the  maximization  [6].  The  algorithm  that  results  is 
very  robust  to  signal  to  noise  ratios  of  2,  and  although  It  Is  more  sensi¬ 
tive  to  model  parameter  choice  It  performs  well  to  signal  to  noise  ratios 
of  1. 

A  recent  paper  by  Nahi  and  Lopez-Mora  [7]  discusses  an  ^d.ternatlve 
and  very  effective  probabilistic  formulation  for  sequential  boundary 
estimation.  In  this  approach  the  likelihood  of  the  picture  data  Is  maxi¬ 
mized  on  an  Individual  row  by  row  basis.  Since  no  attempt  Is  made  to 
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maximize  the  joint  data  likelihood  for  all  of  the  rows  it  is  also  subop- 
timal.  Nor  on  the  surface  does  it  appear  that  there  is  an  easy  way  to 
modify  the  approach  in  order  to  make  optimum  use  of  the  data.  Other 
potential  drawbacks  of  this  approach  are  that  it  does  not  constrain 
boundary  estimates  to  be  continuous  and  hence  it  produces  jagged  boundar¬ 
ies,  and  its  one  dimensional  nature  limits  the  class  of  objects  to  which 
it  is  applicable. 

The  paper  is  organized  as  follows.  In  Section  II  we  present  a  general 
maximum  likelihood  formulation  for  optimal  boundary  estimation.  With  this 
formulation  as  a  guide,  in  Section  III  we  present  a  suboptimal  but  compu¬ 
tationally  tractable  algorithm  for  actual  boundary  estimation.  Two 
alternative  Markovian  boundary  generation  models  are  presented  for  incor¬ 
poration  into  the  algorithm.  Section  IV  demonstrates  the  value  of  the 
likelihood  formulation  for  investigating  algorithm  performance.  The 
suboptimal  algorithm  is  compared  with  both  the  Martelli  algorithm  and 
the  optimal  algorithm  derived  in  Section  II.  We  also  discuss  algoritlim 
performance  in  the  presence  of  artifacts  which  are  inconsistent  with  the 
boundary  model.  Examples  of  eilgorithm  performance  are  given  in  Section  V 
and  some  additional  comments  and  concluding  remarks  are  made  in  Section  VI. 


II.  A  Maximum  Likelihood  Formulation  for 


Boundary  Estimation 

In  this  section  we  present  the  theoretical  framework  on  which  our 
suboptlmal  algorithm  is  based.  In  particular  we  formulate  the  boundary 
estimation  problem  In  terms  of  maximizing  the  joint  log  likelihood  of  an 
hypothesized  object  boundary  and  all  of  the  picture  data  where  the  like¬ 
lihood  Is  derived  by  making  use  of  a  probabilistic  data  generation  model. 

To  begin,  a  digitalized  image  will  be  represented  by  a  picture 
function  gj^j  whose  value  corresponds  to  the  grey  level  of  pixel  {i,j).  We 
model  the  picture  function  g^^  as  consisting  of  two  components — a  true 
picture  component  b^^ ,  and  a  noise  component  n^j ,  so  that  g^^j  =  b^^  +  n^^ . 

The  picture  Is  assumed  to  contain  a  single  object  of  grey  level  r.  lying 

In 

in  a  background  of  grey  level  r  .  where  r.  -  r  ^  =  A  >  0.  For  conveni- 

out  in  out 

ence  we  assume  an  appropriate  constant  has  been  subtracted  from  the  original 

picture  function  so  that  r .  =4  and  r  ^  -  4.  Therefore  b. ,  only  takes 

in  2  out  2  ij  ^ 

on  the  values  y  and  -  y  depending  upon  whether  pixel  (i,j)  lies  inside  or 
outside  the  object.  The  noise  terms  n^^  are  assumed  to  be  Independent, 
identically  distributed  (i.i.d.)  Gaussian  random  variables  with  zero  mean 
and  known  variance,  o^. 

If  an  edge  element  is  defined  as  the  line  segment  separating  two 
adjacent  pixels,  Chen  an  object  boundary  will  consist  of  a  closed  directed 

sequence  of  edge  elements  which  does  not  intersect  itself,  and  will  be 

N  N 

denoted  as  We  make  the  convention  that  a  boundary  Is  gener¬ 

ated  by  moving  in  the  clockwise  direction  around  an  object.  As  depicted 
in  Figures  la  and  lb,  each  directed  boundary  edge  element  1  £  k  £  N, 
can  be  uniquely  described  by  a  threetuple  (l,j,d)  where  1  and  j  correspond 


Co  the  coordinates  of  an  adjacent  pixel,  and  d  to  one  of  the  four  possible 
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edge  directions.  We  can  then  view  the  unknown  edge  sequence  a 

discretely  Indexed  vector  valued  stochastic  process  with  discrete  range 

space.  We  model  this  as  a  Markov  process  of  order  K  with  states 
i+K“l 

and  known  transition  probabilities  Pg(x^|x^  ^).  Appropri¬ 
ate  choices  for  P„(*)  are  given  in  the  section  to  follow. 

Let  us  next  derive  an  expression  for  the  joint  likelihood  L  of  an 
hypothesized  boundary  edge  sequence  and  the  entire  picture  function.  The 
joint  likelihood  L  is  the  product  of  the  likelihood  of  the  hypothesized 

boundary  edge  sequence,  and  the  likelihood  L  t  of  the  picture  function 

P  B 


conditioned  on  the  hypothesized  boundary.  Hence 
S,nL  =  ilnLg  +  ilnLp  j  ^ 

Using  the  Markov  chain  model  we  obtain 


(1) 


JlnLg  =  )lnPg(Xj^)  +  Z  inPg(x^lx^_j^)  +  JlnPj^(N) 


(2) 


where  Pg(Xj^)  is  the  a  priori  probability  of  a  starting  state  Xj^  and  Pj^(N) 
is  the  a  priori  probability  of  a  boundary  being  of  length  N.  Given  a 
specific  boundary  edge  sequence  each  pixel  corresponds  to  either  object 
or  background.  The  Gaussian  noise  model  then  implies  that 


“  ‘■PIB  -  l>^ 

pixel  (i,j) 
in  object 


2,2  'By  2> 

pixel  (i,j) 
in 

background 


(3) 


where  is  a  constant  independent  of  the  choice  of  an  hypothesised 
boundary.  Substituting  (2)  and  (3)  into  (1)  yields  the  following  expres¬ 


sion  for  the  joint  log  likelihood 
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III.  A  Suboptlmal  Boundary  Estimation  Algorithm 

In  this  section  we  present  a  suboptimal  boundary  estimation  algorithm 
which  is  based  upon  maximizing  a  likelihood  function  similar  to  (5),  but 
one  which  can  be  maximized  using  sequential  procedures  such  as  the  A* 
(branch  and  bound)  algorithm  [6].  The  basic  approach  Involves  limiting 
the  use  of  picture  function  data  to  a  small  region  about  an  hypothesized 
boundary. 

As  pointed  out  by  Martelli  [3],  sequential  boundary  estimation  can 
be  viewed  as  finding  a  path  through  a  directed  graph.  Eor  our  formulation 
each  node  of  the  graph  corresponds  to  a  K-dimenslonal  Markov  state  defined 
by  the  boundary  generation  model  introduced  in  the  previous  section.  A 
likelihood  value  is  assigned  to  each  node,  and  it  corresponds  to  the  maxi¬ 
mum  of  the  likelihoods  for  all  paths  leading  to  that  node  from  a  pre¬ 
determined  start  node.  As  shown  in  Fig.  2,  each  node  has  exactly  three 
successor  nodes  corresponding  to  the  three  states  which  can  be  obtained 
by  adding  a  new  edge  element  to  an  hypothesized  boundary  edge  sequence. 

By  defining  a  goal  node  to  be  a  boundary  state  lying  within  a 
neighborhood  of  the  start  node,  as  depicted  in  Fig.  3,  graph  theoretic 
search  procedures  such  as  the  A*  algorithm  can  be  applied  to  obtain  a 
maximum  likelihood  path  from  the  start  node  to  a  goal  node.  The  states 
along  such  a  path  would  then  define  the  boundary  estimate.  These  algo¬ 
rithms  generate  only  the  portions  of  the  graph  needed  to  continue  the 
search.  For  example,  the  basic  step  in  the  procedure  involves  searching 
the  present  graph  for  the  node  with  the  largest  likelihood,  adding  to 
the  graph  any  of  its  three  successors  not  already  on  the  graph,  and 
recalculating  the  likelihood  of  any  of  the  successor  nodes  which  are 
already  on  the  graph.  As  a  result,  it  is  computationally  desirable  to 
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recursively  calculate  the  likelihood  of  a  successor  state  from  knowledge 
of  the  likelihood  of  its  predecessor.  In  view  of  this,  if  is  an 
arbitrary  graph  node  and  j  =  1»2,3  are  its  three  successor  states 

as  depicted  in  Fig.  2,  then  the  following  suboptimal  log  likelihood 
function  can  be  formulated  based  upon  the  optimal  likelihood  (5), 

=  ilnL(x)  +  tnPg(xj^j^|x)  +  j=l,2,3  .  (6) 

D(x^_^l)  is  the  change  in  picture  data  likelihood  caused  by  adding 
node  *£^2^  to  the  boundary  sequence  defined  by  the  most  likely  path  from 
the  start  node  to  node  and  state  transition 

probability  defined  in  the  previous  section.  Methods  for  calculating 
these  two  quantities  are  discussed  in  detail  below. 

Calculating  Picture  Function  Contribution 

Picture  data  is  Incorporated  into  (6)  by  using  pixels  contained  in 
a  swath  about  an  hypothesized  edge  sequence.  Figure  4  gives  examples 
of  data  swaths  for  alternative  boundary  edge  sequences.  This  use  of 
picture  data  is  suboptimal  since  for  the  situations  depicted  in  Figure  4, 
the  log  likelihoods  of  the  two  alternative  paths  would  be  evaluated  using 
different  picture  data.  As  will  be  discussed  in  more  detail  in  the 
section  to  follow,  optimal  use  of  the  data  would  require  comparison  of 
these  two  alternative  boundary  paths  using  likelihoods  which  incorporated 
the  same  picture  data.  Our  data  usage  is  globally  suboptimura,  but  it  is 
locally  optimum.  The  global  suboptimality  is  pertinent  to  the  relative 
frequency  with  which  the  sequential  estimator  leaves  the  vicinity  of  the 
true  path  and  must  backtrack.  The  local  optimality  is  pertinent  to  the 
relative  frequency  with  which  the  estimate  will  be  in  error  by  a  small 
number  of  pixels,  e.g.,  one  or  two  pixels.  Since  the  latter  accuracy 
consideration  may  be  very  important  in  some  applications  we  discuss  it 
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further  below. 

D(x^_l_^)  is  calculated  by  making  use  of  picture  data  contained  in 
a  (4x4)  pixel  block  surrounding  the  last  edge  of  the  boundary  edge 
sequence  leading  to  node  x^.  The  purpose  of  this  (4x4)  template  is  to 
optimally  use  data  in  the  region  of  the  deepest  edge  element  of  a  path 
in  order  to  choose  among  the  following  three  edge  elements  and  to  extend 
the  path  in  an  optimal  way.  This  has  the  advantage  of  permitting  optimal 
detection  for  short  but  arbitrarily  curvy  boundary  subsequences  whereas 
Hueckle  type  operators  [9],  for  example,  are  restricted  to  straight  lines. 

As  shown  in  Fig.  5,  15  possible  extension  edge  sequences  out  of  the 
(4x4)  block  are  considered,  5  for  each  of  the  successor  states 
j=l,2,3.  Assuming  a  boundary  edge  sequence  to  be  generated  by  moving 
in  the  clockwise  direction  about  an  object,  the  original  sequence  and 
any  of  the  15  extension  paths  uniquely  classify  each  of  the  16  pixels 
in  the  block  as  to  being  either  Inside  or  outside  the  object.  Based 
'  upon  (4)  and  (5)  we  have  experimented  with  and  anlyzed  two  alternative 

approaches  for  assigning  likelihoods  to  each  classified  pixel.  In  view 
I  of  (4)  a  natural,  and  our  first  choice  involved  assigning  the  quadratic 

values 

-  I  to(2ira^)  -  -  f)  (7a) 

or 

-  y  iln(2Tra^)  -  -iy(g  +  y)  (7b) 

depending  upon  whether  pixel  (i,j)  was  classified  as  in  object  or  back¬ 
ground  respectively.  However,  as  we  discuss  in  the  Section  IV  to  follow, 
Improved  algorithm  performance  is  obtained  by  assigning  values  more 
analogous  to  the  linear  data  terms  in  (5) .  In  particular  we  assign  the 
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values 

-(A/2a)^  +  (A/2a)(g_/(T)  (8a) 

or 

-{A/2a)^  -  (A/2a) (g^j/a)  (8b) 

depending  upon  whether  pixel  (i»j)  is  classified  as  in  object  or  back¬ 
ground  respectively. 

For  each  node  five  values  k=l,2,...,5  are  calculated; 

one  for  each  of  Lho  five  extension  paths.  Observe  that  those  calculations 
can  be  performed  in  a  computationally  efficient  manner.  It  should  be 
apparent  from  Figure  5  that  for  fixed  J,  each  (4x4)  template,  k,  can  be 


obtained  from  template,  k=l,  by  reclassifying  a  single  pixel.  Hence 
Dj^(x|^l)  can  be  obtained  from  adding  an  appropriate  constant. 

If  any  extension  path  causes  a  loop  with  the  optimal  edge  sequence  leading 
to  state  x^  then  the  corresponding  is  assigned  a  value  of  -«>.  It 

should  also  be  noted  that  some  of  the  data  in  the  (4x4)  template  was  used 


in  calculating  L(x^) .  Hence  only  additions  and  changes  are  incorporated 
into  This  Implies  that  the  complete  data  swath  for  a  boundary 
edge  sequence  consists  of  the  union  of  a  sequence  of  these  (4x4)  blocks. 


Finally  we  define 


D(x^^^l) 


max  \(x-. 
l<k<5  ^  ^ 


Though  the  required  computation  in  calculating  the  15  values  of 
Dk(x^^l)  is  seemingly  large,  it  is  in  fact  modest  since,  as  pointed  out 
above,  these  values  can  be  calculated  recursively.  We  consider  this 
algorithm  to  be  a  very  important  part  of  this  approach  since  the  large 
number  of  template  matchings  required  might  otherwise  be  prohibitive. 
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Calculating  Boundary  Contribution  ?.nPg(x^^j^  |  x^) 

We  have  experimented  with  two  methods  for  calculating  the  Markovian 
transition  probabilities  I  •  To  make  some  initial  comparisons 

with  results  previously  presented  in  the  literature,  we  first  chose 
^B^*l+l^*j^  so  that  -C,nPg(Xj^j^|x^)  roughly  resembled  the  Martelli  curva¬ 
ture  cost  [3].  We  have  also  developed  a  method  for  calculating  transition 
probabilities  based  upon  a  dynamic  boundary  generation  model  and  imple¬ 
mentation  of  a  Kalman  filter. 


The  Martelli  curvature  constraint  [3]  was  based  upon  the  assumption 
that  the  boundaries  of  interest  were  locally  smooth  and  of  low  curvature. 
Under  these  assumptions  we  can  define  our  state  dimension  K=8,  and 

HB^+il^i)=Ae  •’  (10. 

3  f(e  ) 

A  *  T  e  ^  (101 


£(0.)  =  -a|0. I  -  bO:  (10c) 

J  1  J 

where  a>0,  b>0,  are  arbitrary  constants,  and  0.  is  the  angle  depicted 

^  3 

in  Figure  2.  A  is  a  normalization  factor  which  assures  that  7.  Pg(x^^j^  |x^)=l. 

j=l  ^ 

Since  i!.nPg(x^^j^  |x.)  =  <lnA  +  1(0^)  has  a  maximum  for  these  transition 

probabilities  favor  straight  boundaries  and  discourage  sharp  directional 


changes . 

In  implementing  this  scheme  we  observed  that  for  coursely  quantized 
images  Oj  varies  considerably  along  constant  curvature  paths  sucti  as 
circular  arcs.  Hence  to  reduce  tlie  variability  in  our  measurement  of  (L 
we  have  developed  a  look-up  table  procedure  for  determining  directions 
of  line  segments  passing  through  sequences  of  four  edges  such  as  shown 
in  Figure  2.  The  directions  of  the  line  segments  for  the  first  four 


edges  and  for  the  second  four  edges  are  then  compared  to  determine  t'. , 
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The  Kalman  filtering  approach  requires  additional  a  priori  informa¬ 
tion  about  boundary  shape.  It  is  well  known  that  discrete  "time"  linear 
dynamic  systems  driven  by  sequences  of  i.i.d.  random  variables  represent 
Markov  processes.  These  systems  can  provide  good  models  for  many  types 
of  smooth  object  boundaries.  Here  we  will  assume  that  the  boundaries  of 
interest  look  like  perturbed  or  distorted  ellipses.  This  model  is  valid 
for  a  number  of  different  human  organ  boundaries  in  CAT  scans,  or  con¬ 
ventional  tomographs,  and  also  the  cloud  shown  in  Figure  lOe. 

Using  a  Euclidean  coordinate  system,  let  a  sequence  of  points  on 
the  true  object  boundary  be  denoted  by  the  vectors  v(k)  (2x1),  l^k^M. 
Having  restricted  the  image  boundary  to  consist  of  a  sequence  pixel  edges 
t^,  l£i_5N,  points  lying  on  these  edges  represent  quantized  versions  of 
the  points  v(k),  and  will  be  denoted  as  v^(k),  l^k^M.  Since  elliptical 
trajectories  can  be  generated  as  solutions  to  second  order  systems,  an 
appropriate  model  for  distorted  ellipses  would  be 


v(k+l)  =  Av(k)  +  b  + 

(11) 

V  (k)  =  v(k)  +  w.  , 

q  K 

(12) 

where  A  (2x2)  and  b  (2x1)  are  constant  matrices,  and  u  (2x1)  and  w,  (2x1) 

k 

are  sequences  of  zero  mean,  i.i.d.,  Gaussian  random  vectors  with 

2  2 

covariance  matrices  o^I  and  o^I  respectively.  The  process  Uj^  can  be 
viewed  as  the  distortion  mechanism  while  the  process  Wj^  approximates 
the  quantization  error.  We  also  assume  Uj^  and  Wj^  to  be  uncorrelated. 

For  an  ellipse  consisting  of  M  points,  having  a  ratio  of  major  to  minor 
axis  of  p,  a  rotation  angle  of  0,  and  a  center  at  c(2xl),  appropriate 
choices  for  A  and  b  would  be 
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A  = 


Y  +  6( — ^)slnOcosO 

P 

-6  (pcos^6+^in^0) 


6  (-^os^O+psin^O) 

P 

,  2 

Y-6  ( — sinScosO 


(13) 

(14) 

(15) 


b  =  [I-A]c 
2it 

Y  =  cos-^ 

.  .  211 
0  =  sin-^ 

2  2 

From  knowledge  of  A,  b,  o^,  and  0^  one  can  construct  the  well  known 
steady  state  Kalman  filter  (or  predictor)  for  generating  an  estimate 
v(k)  of  v(k) : 

v(k+l)  =  Av(k)  +  b  G[v^(k)  -  v(k) ) 

G  =  AP[a^I+P]"^ 
w 

P  =  APA^  +  0^1  +  APta^l+Pl"^PA^ 
u  w 

To  see  how  this  filter  can  be  used  for  generating  the  transition 

probabilities  Pg(x|^j^Ix^)  consider  Figure  6.  Figure  6  shows  three 
12  3 

possibilities  v  (k) ,  v  (k) ,  and  v  (k)  for  the  point  v  (k) .  These 

q  q  q  q 

12  3 

correspond  to  the  end  points  of  the  three  alternatives,  *'i+l’  *'i+l’ 

for  edge  element  given  t^.  The  estimate  v(k)  of  v(k)  is  also  shown. 

Having  assumed  Uj^  and  Wj^  to  be  Gaussian,  one  can  show  that  if  v^(k)  were 

unconstrained  then  it  would  be  Gaussian,  and  conditioned  on  v^(k-l)  it 
.  2 

would  have  mean  v(k)  and  covariance  P  +  a  I.  However,  since  v  (k)  is 

W  Q 

constrained  to  lie  on  a  pixel  edge  we  can  discretize  the  Gaussian  density 
function  to  obtain  approximate  transition  probabilities. 

Let  the  Markov  process  state  x^  be  the  single  edge  element  t^ , 


P(v-^(k+l)|v  (k)) 


i  P„(vJ(k),v(k),afl+P)/  7.  P_(vj(k),v(k),ajl+P) 

O  Q  W  .  ,  u  Cl  W 

s  J  =  1 


(16) 


then 
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where 

Pp(x,m,>:)  =  (2it(>;|)  ^exp(-  y(x-m)^>:~^(x-in))  (17) 

In  Implementing  this  scheme,  since  the  dynamic  system  does  not 
generate  equidistant  point,  we  have  found  it  helpful  to  choose  M=3N  so 
that  on  the  average  we  generate  3  estimates  per  pixel.  To  obtain  the 
value  v(k)  used  in  (17)  the  predictor  is  run  sequentially  from  v(k-l) 
until  a  pixel  boundary  is  crossed. 

Graph  Searching 

The  state  likelihoods  calculated  by  (6)  are  used  in  conjunction 
with  a  modified  A*  algorltliin  to  find  the  boundary  estimate.  First 
since  our  (4x4)  data  template  was  developed  under  the  hypothesis  that 
the  boundary  passed  through  the  template,  before  we  expand  a  state,  i.e. 
add  its  successor  states  to  the  graph  we  test  this  hypothesis.  In 
particular,  if  the  probability  that  the  (4x4)  template  (centered  about 
the  last  edge  element  of  the  state  to  be  expanded)  is  completely  Inside 
the  object  or  background  is  above  a  fixed  threshold  then  this  state  is 
not  expanded.  Second,  to  make  the  algorithm  computationally  feasible  in 
terms  of  both  computer  memory  and  execution  time  requirements  it  is 
necessary  to  periodically  prune  the  graph  of  nodes  which  had  little 
probability  of  being  on  a  boundary.  If  addition  of  a  set  of  successor 
states  j=lf2,3  increased  the  maximum  depth  of  the  graph,  then  all 

graph  branches  a  distance  (T+1)  or  more  back  from  the  states 
are  pruned  from  the  graph  provided  they  do  not  lie  on  the  most  likely 
path  from  the  starting  state  to  states 

In  Section  V  to  follow  we  present  some  examples  to  illustrate  the 


algorithms  performance. 
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IV.  Interpretation  and  Analysis 

In  this  section  we  show  how  the  mathematical  framework  introduced 
in  Section  II  can  be  used  to  do  some  simple  but  Insightful  analyses  of 
the  boundary  estimation  algorithm  outlined  in  Section  III.  The  first 
analysis  compares  the  operation  of  this  algorithm  with  that  of  the 
Martelli  algorithm,  while  the  second  examines  the  effects  of  the  algo¬ 
rithms  restricted  use  of  picture  function  data  and  shows  that  performance 
will  not  be  significantly  impaired.  Finally,  the  third  analysis  shows 
how  it  is  possible  to  predict  algorithm  performance  when  certain  unexpected 
artifacts  are  encountered. 

Comparison  with  Martelli  Algorithm 

If  the  first  method  for  calculating  transition  probabilities,  using 
local  curvature  Information,  is  incorporated  into  the  algorithm,  the 
essential  difference  between  it  and  the  Martelli  algorithm  is  the  method 
in  which  picture  function  data  is  employed  for  estimation.  The  Martelli 
algorithm  searches  a  directed  graph  for  a  minimum  cost  path  where  one 
component  of  the  nodal  cost  function  is  a  directional  derivative  calcu¬ 
lated  from  local  picture  function  data.  This  is  somewhat  analogous  to 
the  picture  function  component  of  (6)  obtained  by  classifying  data  in 
the  (4x4)  template.  In  this  subsection  we  compare  the  performance  of  the 
two  algorithms  by  looking  at  pairs  of  edge  sequences  (paths) .  One  path 
is  assumed  to  lie  on  the  true  boundary  while  the  other  is  assumed  to  lie 
in  the  object  or  background.  We  then  calculate  an  error  probability  for 
each  algorithm, that  is  the  probability  that  the  path  lying  in  the  object 
or  background  is  more  likely  or  of  least  cost.  We  first  consider  the 
case  of  a  single  erroneous  edge  and  show  that  our  new  algorithm  has  a 
lower  probability  of  initially  leaving  a  correct  path.  The  implication 
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of  this  Is  that  fewer  nodes  will  be  added  to  the  graph.  We  next  look 
at  longer  paths  where  one  path  is  straight  and  the  other  curvy.  Here 
it  is  shown  that  our  algorithm  performs  better  when  the  true  boundary  is 
curvy  while  Martelli's  algorithm  performs  better  when  the  true  boundary 
is  straight.  In  fact,  in  the  latter  case,  the  Martelli  algorithm  per¬ 
forms  better  than  the  optimal  algorithm  which  maximizes  (5).  This  can 
be  explained,  however,  by  realizing  that  the  optimal  formulation  assumes 
no  a  priori  information  about  boundary  curvature  while  the  Martelli 
algorithm  is  designed  to  favor  straight  paths. 

First  consider  the  patch  of  picture  in  Figure  7.  It  shows  a  correct 
edge  sequence  (lying  on  the  object  boundary),  along  with  an  erron¬ 

eous  edge  ty  Let  us  calculate  the  probability  of  each  algorithm  expand¬ 
ing  (generating  successor  states  to)  the  incorrect  state  x '={..., tj^,t2, 
tj}  rather  than  the  correct  state  x={ . . . ,tj^,t2»t2^*  Assuming  the  decision 
is  based  solely  on  use  of  picture  function  data,  Martelli's  algorithm 
would  decide  x'  rather  than  x  if  the  data  cost  c(x')  were  less  than  the 
cost  c(x).  The  cost  function  used  by  Martelli  approximates  the  direc¬ 
tional  derivative  across  an  edge  and  makes  use  of  4  pixels  in  the  row 
(or  column)  containing  the  edge  (2  on  each  side) .  In  addition,  the 
constant  2A  is  added  to  Insure  that  the  cost  is  zero  if  calculated  across 
a  true  boundary  in  the  noise  free  '-ase.  Using  this  function  we  obtain 

c(x  )  =  2A  +  g2j  +  g22  -  ^23  ~  ^24 

c(x)  =  2A  +  gj^2  +  g2j  -  gj^  -  g^j 

Therefore 

P(deciding  x'|x  correct) 

=  P(c(x)-c(x')  >  0|x  correct) 

"  P(2g23  +  ^13  ®24  ~  ^21  ”  ^22  ~  ^33  ~  ^43  ■  Oix  correct 

=  P(G  >  0|x  correct) 


Conditioned  on  x  being  correct*  G  is  a  normally  distributed  random 

2  2 

variable  with  mean  -A  and  variance  lOo  ,  i.e.,  G''N(-A,10a  ),  This  implies 


P(G>0|x  correct)  =  /  (2n)  exp(-s  /2)ds 


/TcT 


Using  our  algorithm  calculation  of  the  data  likelihood  for  both  x 
and  x'  involves  classifying  the  same  16  pixels  shown  in  Figure  7.  The 

only  difference  between  the  two  likelihoods  is  the  value  assigned  pixel 

'  2 

(2,3)  for  the  likelihood  ilnL(x)  it  is  assigned  -(A/2a)  -  (A/2(>)  (823/^) 

2 

while  for  !lnL(x')  it  is  assigned  -(A/2o)  +  ( A/2a)  (g23/o)  .  Therefore 

P(deciding  x'|x  correct)  =  P( !lnL(x'HnL(x)  >  0|x  correct) 

2  I 

=  P((A/o  )823  ^  correct) 

=  P(H  >  0[x  correct) 

2  2 

Conditioned  on  x  being  correct  H  is  N(-(l/2) (A/a)  ,(A/o)  ) 


P(H>0|x  correct)  =  /  (2ti)  ^exp(-s  /2)ds 


Since  the  lower  limit  of  the  integral  (19)  is  1.6  times  that  of  (18)  the 
latter  probability  is  higlier.  As  a  numerLcal  example,  consider  the  case  of  a 
signal  to  noise  ratio,  (A/a),  of  one.  Then  the  lower  integral  limits  in  (18) 
and  (19)  are  .5  and  .625,  and  the  corresponding  error  probabilities  are  .265 
and  .309  respectively.  This  gives  an  Indication  that  tliis  algoritlim  is  less 
likely  to  explore  extraneous  paths  than  Martelli's,  or  equivalently,  it  would 

return  substantially  the  same  performance  in  the  presence  or  larger  noise 
variance. 

Next  consider  the  eight  pairs  of  paths  shown  in  Figure  8.  In  each 
case  path  1  is  straight  and  path  2  is  a  curvy  perturbation  obtained  by 
moving  four  pixels  from  inside  object  to  inside  background  or  vice-versa. 

If  paths  1  and  2  have  Martelli  costs  c^  and  and  log  likelihoods  ^nL^^ 
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and  ilnL2  defined  by  our  algorithm,  then  the  algorithms  choose  paths  1 
or  2  in  accordance  with 
path  2 

G  =  (c.-c_)  0  (20a) 

path  1 

path  2 

H  =  (J?.nL--!lnL,)  ^  0  .  (20b) 

path  1 

G  and  H  will  each  have  two  components  a  boundary  statistic  denoted  as 
B, either  curvature  cost  or  transition  probability,  and  a  data  statistic 
denoted  as  A,  either  gradient  or  template.  Conditioned  on  one  path  being 
correct  (lying  on  boundary)  and  the  other  being  incorrect,  then  only  A  is 
a  random  variable  and  the  total  statistic  (A+B)  is  normally  distributed 
with  mean  and  variance  which  are  a  function  of  which  path  is  assumed 
correct.  Hence  the  probability  of  an  error,  that  is  each  algorithm 
finding  the  incorrect  path  of  least  cost  or  more  likely,  is  of  the  form 
Pp  =  P(A+B  ^  0 1  paths  1  and  2) 

00  J 

=  J  (2ir)  ^exp (-s^/2)ds  (21) 

P 

where 

D  =  _ +(B+Mean(A)) _  ,22'> 

Standard  deviation(A)  ’  ^  ■' 

Assuming  B=0  and  a  signal  to  noise  ratio  of  one.  Table  1  gives  p  and 

for  each  of  the  eight  cases  both  when  the  straight  path  1  is  assumed 

correct  and  also  when  the  curvy  path  2  is  assumed  correct.  When  both 

possibilities  are  equally  likely  tlien  tiie  average  probability  of  an 

error,  P  ,  is  just  the  average  of  P  for  the  two  cases.  This  is  also 
a  p 

calculated  in  Table  1.  If  an  algorithm  were  employed  which  chose  paths 
according  to  the  optimal  likelihood  function  (5)  then  the  only  pixels 
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entering  the  error  statistic  would  be  the  four  lying  between  the  two  paths. 
Hence  one  can  easily  calculate  p  and  for  this  algorithm  and  this  data 
is  given  in  Table  1  as  well.  A  common  property  of  our  algorithm  and  the 
optimal  is  the  symmetry  of  the  error  statistic.  One  obtains  the  same 
error  probability  regardless  of  whether  path  1  or  path  2  is  assumed 
correct.  This  is  not  true  for  the  Martelli  algorithm.  In  particular 
the  Martelli  algorithm  does  considerably  better  when  the  correct  path  is 
straight.  In  fact  in  this  case  it  out  performs  the  optimal.  This  is 
because  it  has  been  designed  to  favor  straight  paths,  while  neither  the 
data  statistic  for  our  algorithm  nor  the  optimal,  favors  any  type  of 
curvature.  On  the  other  hand,  because  of  the  assymetry  of  the  Martelli 
statistic  our  algorithm  gives  considerably  better  performance  if  the 
correct  path  is  curvy.  Although,  as  expected,  it  does  not  perform  as 
well  as  the  optimal.  A  more  detailed  comparison  of  our  algorithm  and 
the  optimal  is  given  in  the  next  subsection.  Finally  in  comparing  overall 
performance,  observe  that  for  five  of  the  cases  our  algorithm  gives  sig¬ 
nificantly  lower  average  error  probability  while  Martelli ’s  algorithm 
does  significantly  better  iit  only  two  of  the  cases. 

Justification  for  the  Suboptimal  Boundary  Estimator 

In  Section  IT  we  derived  a  likelihood  function  (5)  for  optimal  MAP 
estimation  of  the  state  sequence  Ix^  I  or  equivalently  the  boundary  edge 
sequence  {t^}.  A  key  feature  of  (5)  is  the  way  it  incorporates  picture 
function  data.  First,  it  makes  use  of  all  the  picture  function  data. 
Second,  as  seen  in  (A),  the  data  enters  quadrat ically.  However,  after 
simplification  only  the  linear  terms  effect  the  maximization.  To  develop 
a  computationally  feasible  algorithm  picture  data  had  to  be  used  in  a 
suboptimal  manner.  In  this  subsection  we  Investigate  analytically,  per¬ 
formance  effects  caused  by  our  suboptiin.il  data  usage.  In  particular  wc 
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show  that  our  boundary  estimator  performs  substantially  the  same  as  the 
optimal.  We  also  give  justification  for  using  picture  data  linearly  in 
our  suboptimal  algorithm  rather  than  quadrat ically. 

As  done  in  the  previous  subsection,  we  consider  two  alternative  paths. 
Path  1  is  assumed  to  lie  on  the  true  boundary  and  we  calculate  the  proba¬ 
bility  of  an  error,  that  is  that  path  2  appears  more  likely.  Rather  than 
consider  specific  paths  such  as  in  Figure  8,  we  will  consider  two  more 
general  cases.  As  shown  In  Figure  4a,  in  case  1  we  will  assume  that  the 
two  paths  are  far  enougli  apart  so  that  tlic  two  data  swaths  used  by  the  sub- 
optimal  algorithm  do  not  Intersect.  In  case  2,  Figure  4b,  we  assume  that 
the  paths  are  close  enough  together  so  that  the  swath  sections  to  the 
left  of  path  2  and  to  the  right  of  path  1  overlap. 

For  each  case  we  consider  three  algorithms.  The  optimal  algorithm 
based  on  (5)  makes  use  of  all  the  data  in  the  picture  region  containing 
the  two  paths.  However,  a  more  practical  approach,  although  less  optimal, 
is  to  incorporate  only  the  data  in  the  two  swaths,  but  use  all  of  it  in 
calculating  the  likelihood  for  each  path.  We  compare  this  near  optimal 
algorithm  with  our  algorithm  and  it  will  be  designated  Algorithm  1. 

The  second  algoritlira  considered  is  the  final  version  of  our  suboptimal 
algorithm  where  picture  data  enters  the  likelihood  function  (6)  linearly 
according  to  (8).  With  this  algorithm  the  likelihood  of  a  specific  path 
contains  only  data  from  the  swath  about  that  path.  As  mentioned  in 
Section  III,  since  picture  data  enters  equation  (4)  quadratically,  in 
our  early  experimental  work  we  chose  to  incorporate  picture  data  into  (6) 
quadratically  according  to  (7).  In  order  to  demonstrate  the  improved 
performance  obtained  by  using  the  linear  rather  than  the  quadratic  terms, 
Algorithm  3  is  also  analyzed.  It  is  the  same  as  Mgorithm  2  but  incor¬ 
porates  picture  data  using  (7)  rather  than  (8). 
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Before  proceeding  some  additional  notation  is  needed.  Consider 
Figure  4  and  let  i!.e(1,2),  m£(L,R)  denote  the  pixel  subset  con¬ 

tained  on  side  m  of  swath  J.  (swath  containing  path  1)  .  Next  let 
denote  the  entire  set  of  pixels  in  swath  9.  so  that  S  =  S  U  S  .  Define 

)i  Jclj  2,R 

the  number  of  pixels  in  set  or  as  N^^^  or  respectively.  Finally 

for  the  overlapping  case  of  Figure  4b,  let  S  =  S,  S-  and  N.  be  the 

1  XK  L 

subset  and  number  of  pixels  in  the  intersection  of  S,  and  S„  . 

XR  ZL 

Case  1  -  Algorithm  1:  In  this  case  we  assume  and  to  be  dis¬ 
joint,  and  path  1  lies  on  the  true  boundary.  Algorithm  1  uses  likelihood 
(5)  but  only  for  pixel  data  in  and  $2.  Defining  and  L2  as  the  like¬ 
lihoods  of  paths  1  and  2  respectively,  then  an  error  is  made  if  >  0 
where 


Qfl  =  {.nL2  -  f-nLj^ 


r.  in  P  (x Jx  ,)  -  r.  in  P  (x  lx  ) 
path  2  path  1 


(l,j)cS 


(A/2o) (g. ./o) 


>:  (A/2o)(g,,/o)  (23) 


2L 


(i,j)cS 


IR 


As  in  the  previous  subsection  =  A^^^  +  B  where  B  is  the  boundary 
statistic  consisting  of  the  transition  probability  terms  in  (23),  and 
Aj^j^  the  data  statistic  consisting  of  the  picture  function  terms.  Con¬ 
ditioned  on  two  specific  paths  only  A^^^^  is  a  random  variable.  It  is 
Gaussian  with  distribution 

A^^  '  N(-2(Nj^j^+N2L)(A/2a)^,  4(N^j^+N2j^)  (A/2o)^)  • 

Hence,  conditioned  on  paths  1  and  2  the  probability  of  error 
P^j^=P(Qj^j^  >  0 1  paths  1  and  2)  = 

00 

“  /  (2tt)  '^exp(-s^/2)ds  (25) 

•’ll 

where 

-B-Mean(Ajj^) 

*^11  Standard  deviation (A^j^) 


(26) 
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For  B=0  and  we  obtain 


'll 


/2 


(27) 


Numerical  values  for  and  when  the  signal  to  noise  ratio 


(A/o)  is  one  and  two,  and  when  the  number  of  pixels,  N,  in  swaths  S 
and  is  4  and  8,  are  given  in  Table  2. 


2L 


Case  1  -  Algorithm  2:  For  our  suboptimal  algorithm  the  likelihood 
associated  with  path  1  only  contains  picture  data  in  swath  while  tliat 
associated  with  path  2  only  contains  data  from  swath  2.  Hence  using 


equations  (6)  and  (8) ,  if  and  are  the  likelihoods  associated  with 
paths  1  and  2  respectively  then  an  error  is  made  if  ^  0  where 


Q^2  “  ~  ~ 


^  JnPg(x  |x._^)  -  j:  enPgCx  lx 
path  2  path  1 


1  (A/2o)(g  /o)  + 

(i,j)£S^^  (i,j)  S 


Z  (A/2a)(g^j/o) 
2R 


Z  (A/2o)(g,,/o)  -  Z  (A/2o)(g.,/o) 


2L 


+  (N^-N2)(A/20)‘ 


(28) 


^12  ~  ^12  ®  '^here  B  is  the  boundary  statistic  and  Aj^2  the 


data  statistic,  then 

Aj^2  '  N(-2N2L(A/2a)^  ,  (N^  + 

Again  we  can  calculate  an  error  probability 
^12  ^  ^^*^12  ^  0 1  paths  1  and  2) 


(29) 


=  J  (2tt)  ^exp(-s^/2)dE 
"12 


(30) 


For  B=0  and  =  N2j^  =  N2j^  =  N  we  obtain 

»12  ■ 

Numerical  values  for  Pj^2  and  Pj2  can  be  found  in  Table  2.  Observe 
that  ^12’  t»ur  suboptimal  data  usage  results  in  a 


(31) 
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boundary  estimator  which  is  quite  good.  In  fact  for  as  few  as  8  pixels 
in  a  swath  (N=4) ,  and  a  signal  to  noise  ratio  (A/a)  =  1,  the  error  proba¬ 
bility  is  only  .150.  Also,  as  is  apparent  from  examination  of  Table  2, 
if  we  consider  longer  paths  with  more  pixels  in  a  swath,  or  pictures 
with  larger  signal  to  noise  ratios,  the  error  probability  quickly  decreases 
to  zero. 


Case  1  -  Algorithm  3:  If  we  use  an  algorithm  similar  to  algorithm 
2  but  one  where  data  likelihoods  are  calculated  using  the  quadratic 
terms  (7)  rather  than  the  linear  terms  (8)  we  obtain  an  error  statistic 
Qi3  =  A33  +  B  where 


Ai3  - 


r-  1  ^  ^  .  ,,  1  .  h.2 

„  2  ^®ij  2^  2  ^®ij  "  2^ 


(i,j)cS2j^  2o 


1  .  .  A.2  _  1  ,  A.2 

2  ^®ii  2^  ~  ^  2  ^®ii  ~  2 


+  j  (Nj^-N2)lln(2Tta^) 


variance 


Mean(A^3)  =  |(Nj^-N2)  (l+Pn(2T..T^))  - 
=  ■|(N^+N2)  +  4N2j^(A/2o)^ 


(32) 


Although  A33  is  not  Gaussian  if  and  N2  are  large  then  it  is 
reasonable  to  approximate  it  as  Gaussian.  Since  4^3  lias  mean  and 


For  B=0  and  N  as  defined  above 


^  ^  (H-2(^)^)  -' 


(33) 


Numerical  values  for  and  the  error  probability  P33  under  the 


Gaussian  assumption  are  given  in  Table  2.  Clearly  Algorithm  2  outperforms 
Algorithm  3  for  all  finite  signal  to  noise  ratios,  and  significant  improve¬ 
ment  is  obtained  for  signal  to  nolsi*  r.it  ios  of  oiu'  or  loss. 
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Case  2  -  Alt;orithin  1;  As  illustrated  in  Figure  (4b)  for  this  case 
we  assume  paths  1  and  2  are  closer  together  so  that 

empty.  However  we  do  not  allow  S,  to  intersect  S-  or  S.  to  intersect 

X  ZR  2. 

The  overlap  docs  not  influence  the  performance  of  Algorithm  1 
hence  and  the  error  probability 

Case  2  -  Algorithm  2:  Although  for  algorithm  2  the  error  statistic, 
Q22»  is  still  given  as  in  Case  1  by  equation  (28),  the  overlapping  data 
in  and  82^^  effect  the  distribution  of  the  corresponding  data  statistic 
A22.  The  pixels  contained  in  the  overlapping  set  Sj  enter  twice  rather 
than  once  causing  to  have  a  larger  variance  than  A^2-  i*'  particular 

A22  --  N(-2N2l(A/2o)^,  (Nj^+N2+2Np(A/2o)^)  (34) 

This  causes  a  degradation  in  performance  in  comparison  with  both  Case  I, 
Algorithm  2,  and  Cases  1  and  2,  Algorithm  L.  For  example,  if  we  consider 
a  limiting  situation  where  =  S2j^  so  that  =  N  we  obtain 


P22  «)■<?) 


Hence  0^2  ~  (^/3) ^ 


^22*  ^11 


P22  *^3  f>22- 


(35) 

However,  as  illustrated 


in  Table  2,  performance  is  still  quite  good. 

Case  2  -  Algorithm  3;  For  Algorithm  3  the  data  overlap  actually 
causes  an  improvement  in  performance.  Tn  this  case  the  non-Gaussian 
data  statistic  A2J  has 


Mean  (k^^)  =  |■(NJ^-N2)  (l+((.n(2«o^) )  -2N2j^(A/2n)^ 

Var  (A23)  =  |(N^+N2-Nj)  +  4N2l(A/2o)^ 

Making  the  Gaussian  approximation  and  taking  =  N  we  obtain 


(36) 

(37) 


P  =  |(N)^2(^) 

^  (l+^(^)‘) 
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3/CL\2n?5 


2'a' 


(38) 


25 


Although  the  overlap  causes  a  degradation  in  the  performance  of  algorithm  2 
and  an  improvement  in  that  of  Algorithm  3,  as  can  be  seen  from  Table  2, 
for  the  case  of  small  signal  to  noise  ratios.  Algorithm  2  still  outper¬ 
forms  Algorithm  3.  This  is  the  more  important  situation  since  for  large 
signal  to  noise  ratios  all  error  probabilities  are  small. 

Estimating  Projections 

We  would  like  to  be  able  to  predict  algorithm  performance  when  the 
object  boundary  contains  artifacts  which  are  somewhat  Inconsistent  with 
the  boundary  generation  model.  Consider,  for  example  the  projection  shown 
in  Figure  8,  and  assume  the  notation  and  terminology  of  the  previous 
analyses.  If  the  state  transition  probabilities  Pg(’l’)  small  for 

paths  of  high  curvature,  that  is  the  Markov  process  boundary  model  favors 
smooth  straight  paths,  then  the  contribution  to  the  boundary  statistic  B 
for  path  1  will  be  larger  than  that  for  path  2.  Thus  the  statistic  B 
will  be  positive,  and  tend  to  support  an  erroneous  decision  in  favor  of 
path  2.  The  only  way  that  the  overall  decision  statistic,  A+B,  can  support 
a  correct  decision  is  for  A+B<<0,  or  A  sufficiently  negative.  This  will  be 
the  case  if  and  N2  are  large  enough,  that  is  if  there  is  sufficient 
data  so  that  its  contribution  to  the  decision  statistic  outweighs  the 
boundary  model  contribution.  More  specifically  we  require 

p  =  — -B  -Mean(A) - 

Standard  deviation(A)  ^  ' 

Using  (29)  we  obtain 

-B  +  2N  (A/2o)^ 

- £li - >>  1  (40) 

(N^+N2)  (A/2a) 

If  the  projection  is  long  relative  to  the  number  of  edge  elements  used 
in  calculating  transition  probabilities  for  B,  B  will  remain  roughly 
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constant  as  and  N2  increase.  Using  (40),  and  estimating  B  for  a 
specific  boundary  model  allows  one  to  determine  roughly  the  required 
swath  width  for  correct  estimation,  A  more  accurate  determination  of 
the  estimators  ability  to  track  around  projections  requires  a  more 
careful  analysis  of  the  sequential  behavior  of  the  algorithm.  An  example 
illustrating  the  results  of  applying  our  algorithm  to  a  projection  such 
as  in  Figure  9  is  given  in  the  following  section. 

An  alternative  approach  to  handling  boundaries  with  unusual  artifacts 
such  as  occasional  projections  or  others,  is  to  treat  these  as  pattern 
classes  to  be  recognized. 


L 


V .  Examples 


In  this  section  we  present  some  examples  of  the  boundary  estimates 
computed  by  our  algorithm.  The  images  we  consider  are  shown  in  Figures 
lOa-e.  Figures  lOa-c  are  artificially  generated  noisy  images,  while 
Figure  lOd  is  a  FLIR  image  of  a  tank  and  Figure  lOe  is  a  satellite  image 
of  a  cloud. 

First  consider  Figure  10a.  It  shows  a  perfect  ellipse  in  an  additive 
Gaussian  noise  field  resulting  in  a  signal  to  noise  ratio  A/a  =  1.  Fig¬ 
ure  11a  shows  both  the  actual  elliptical  boundary  and  tlie  estimated 
boundary  when  transition  probabilities  were  calculated  using  the  local 
curvature  algorithm.  The  parameters  a  and  b  in  equation  (10)  were  2.0 
and  5.0  respectively.  For  the  estimate  shown  in  Figure  11b  the  transi¬ 
tion  probabilities  were  calculated  by  use  of  the  Kalman  filtering  algorithm 


0.568 

•0.0365 


0.0138 

0.5683 


As  expected,  since  the  Kalman  filter  approach  makes  use  of  more  ■: 


information,  it  generates  better  estimates.  Computationally,  it  is  also 
more  efficient  in  terms  of  the  number  of  nodes  placed  on  the  graph.  The 
elliptical  boundary  contains  roughly  180  edge  elements.  Hence,  if  an 


algorithm  did  not  expand  any  extraneous  graph  nodes  it  would  expand 
slightly  fewer  than  180  nodes.  The  curvature  algorithm  expanded  307 
nodes  while  tlx*  K.-iIm.m  filter  algorithm  expaiHle<t  only  200  nodes. 
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Next  consider  Figure  10b.  It  shows  a  distorted  ellipse  whose 
boundary  was  generated  by  use  of  the  dynamical  model  given  In  (11)  and 
(12) .  It  also  Is  Imbedded  In  an  additive  Gaussian  noise  such  that 
A/o  =  1.  Figure  12a  shows  the  true  and  estimated  boundary  using  the  Kalman 
filter  algorithm,  and  the  state  space  model  parameters  which  generated 
the  true  boundary.  For  this  model 


p.o 

0.0209 

[-0.058 

1.0 

p0.6688“ 

LI.856 

0.9285 

0.0194 

^-0.0538 

0.928^ 

These  matrices  were  generated  by  choosing  N=180,  M=3N,  0=90®, 

T 

p=3/5,  c  =  [32,32],  o  =1  and  a  =l/t/i^.  To  demonstrate  the  robustness 

u  w 

of  the  algorithm  to  poorly  estimated  model  parameters,  we  perturbed 

these  parameters  by  roughly  20%.  Specifically  we  chose  N=160,  0=70°, 

p=3/A,  c"^  =  [26,26],  o  =1,  and  a  =l//i2.  This  results  In 

u  w 

r 1.0  0.031?] 


A  = 


[-0.0489  1.0 


b  = 


-0.962 

1.2714 


K  = 


""0.9285 

-0.0454 


Figure  12b  shows  the  estimated  and  true  boundaries  In  this  case.  While 
the  resulting  performance  In  the  two  cases  Is  very  similar.  It  appears 
that  accuracy  of  model  parameter  estimation  can  substantially  influence 
computation  time.  With  the  correct  model  parameters  the  algorithm 
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expanded  only  190  nodes  while  with  the  incorrect  parameters  the  algorithm 
expanded  438  nodes. 

Figure  10c  shows  a  circular  object  with  a  rectangular  projection  of 
the  type  discussed  in  Section  III.  The  signal  to  noise  ratio  in  this 
picture  is  also  one.  Figure  13  shows  the  true  and  estimated  boundary  for 
this  picture  using  the  curvature  algorithm  for  calculating  transition 
probabilities.  In  this  case,  the  data  swaths  used  in  the  likelihood 
calculations  were  large  enough  to  assure  proper  performance. 

Finally,  Figures  14  and  15  show  the  results  of  applying  our  algorithm 
to  real  image  data.  Figure  14  shows  the  estimated  boundary  for  the  FLIR 
image  of  the  tank  shown  in  Figure  lOd  while  Figure  15  shows  the  estimated 
boundary  for  the  satellite  cloud  image  of  Figure  lOe. 
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VI.  Comments  and  Conclusions 


By  viewing  boundary  estimation  in  terms  of  MAP  estimation  of  a  Markov 
process  state  sequence,  in  Section  II  we  presented  a  maximum  likelihood 
framework  to  serve  as  a  basis  for  the  design  and  analysis  of  sequential 
boundary  estimation  algorithms.  While  optimal  estimation  required  use  of 
all  the  picture  data,  computational  constraints  limited  the  amount  of 
data  which  could  be  practically  employed  by  an  algorithm.  As  a  result, 
the  particular  subopt Imal  algorithm  presented  in  Section  III  only  incorpor¬ 
ated  picture  data  in  swaths  about  hypothesized  boundary  edge  sequences. 

The  results  of  Section  IV  demonstrated  the  power  of  the  likelihood  formula¬ 
tion  in  doing  formal  algorithm  analysis  for  the  purposes  of  comparing 
different  algorithms  and  improving  performance  of  a  specific  algorithms. 

It  should  be  pointed  out  that  the  particular  algorithm  presented  in 
Section  III  is  but  one  of  many  possible  suboptimal  algorithms  which  are 
consistent  with  the  general  maximum  likelihood  formulation,  and  there  is 
still  potential  for  developing  algoirthms  with  both  improved  performance 
and  reduced  computational  complexity.  Design  of  such  algorithms  involves 
two  important  steps.  One  must  choose  an  appropriate  Markovian  boundary 
generation  model,  and  one  must  decide  on  how  to  incorporate  picture  data 
into  the  likelihood  computations.  With  this  in  mind  we  make  some  final 
comments  on  these  subjects. 

Boundary  Model  Design 

In  choosing  a  boundary  model  one  must  decide  upon  the  type  of  a  '.n'lori 
information  regarding  boundary  curvature  and  shape  which  is  available. 

This  information  can  be  generally  classified  as  local  or  global  in  nature. 
For  example  smoothness  is  a  local  property  of  a  boundary,  while  closure 
is  a  global  property.  Though  the  two  models  we've  made  use  of  incorporate 
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both  local  and  global  information,  the  curvature  model  emphasizes  local 
information  while  the  dynamical  system  model  emphasizes  global  information. 
In  particular,  the  curvature  model  is  described  in  terms  of  behavior  of 
short  sequences  of  edge  elements.  It  does  not  make  use  of  the  fact  that 
the  boundaries  of  interest  are  closed.  On  the  other  hand,  although  curva¬ 
ture  information  is  implicitly  incorporated  into  the  parameters  of  the 
dynamic  model,  this  model  was  explicitly  developed  to  exploit  a  prior'i 
knowledge  that  the  boundaries  of  interest  were  closed. 

Picture  Data  Usage 

In  view  of  the  optimal  formulation  one  is  clearly  interested  in  making 
use  of  as  much  picture  data  as  possible.  However,  to  take  advantage  of 
the  computationally  attractive  sequential  maximization  procedures  such  as 
A*  or  dynamic  programming,  one  must  limit  the  use  of  picture  data.  We 
chose  to  consider  (4x4)  data  windows  and  hence  limited  our  use  of  picture 
data  to  swaths  of  roughly  2  pixel  width  on  each  side  of  an  hypothesized 
boundary.  Although  as  demonstrated  in  the  comparison  with  the  Martelli 
algorithm  in  Section  IV,  this  window  is  useful  in  preventing  a  first 
erroneous  edge  element  decision,  since  different  swaths  contain  different 
data,  performance  is  degraded  In  comparison  with  optimal  data  usage.  As 
a  result,  a  computationally  feasible  algorithm  which  makes  use  of  a  larger 
data  window  containing  a  greater  number  of  paths  for  comparison,  could 
considerably  improve  performance.  Use  of  a  larger  window  also  appears 
Important  for  Images  containing  large  numbers  of  pixels  where  grey  level 
changes  along  boundaries  may  be  more  gradual.  One  approach  for  sequentially 
estimating  boundaries  through  larger  window  blocks  is  given  in  [5].  It 
might  also  be  mentioned  that  a  more  detailed  discussion  of  the  dependence 
of  boundary  error  on  data  usage  is  given  in  [1]. 
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In  conclusion,  the  guiding  principle  in  this  work  has  been  to  formu¬ 
late  an  optimal  approach  to  the  problem,  and  then  to  design  a  subopt imal 
realization  whose  performance  can  be  fairly  well  understood  and  analyzed. 
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Martelli  Algorithm 


Table  1  -  Comparison  of  Martelli,  Daboptimal,  and  Optimal  Algorithms 


Figure  ib.  lixample  of  a  closed 
directed  boundary. 


Figure  2.  Relation  between  boundary  edg.e  sequences 
and  graph  nodus  (Markov  states). 


boundary 


Figure  3.  Graph  structure  corresponding  to  two 
alternative  object  boundaries. 
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Fig,  /4a  Examples  of  alternative  boundary  path  segments 
with  disjoint  data  swaths. 
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Fig.  4b  Example  of  paths  with  overlapping  data  swaths. 


paths  for  calculating  data 


Kalman  filter  estimate  v(k)  and 
corresponding  quantized  data  points 

Vq(k). 


(4xA)  section  of  a  picture. 


^tundiiS  atiiua’* 


45 


oaooooo«»««seae*c 


»  c«  « 

*  k  3 

*30 

s :] 

3  «• 

3  * 


t  men  ******  * 


«  a«3 

*  ft  *  » 

oo  * 

3C 


a  * 

3  * 

O  ♦ 

D#- 

a 

a 

a 

3  * 
ft 
3  * 
3  # 


a 

♦  3 
C  * 

3*  * 
30* 

C« 


a* 

aa«#  _ 
oasa 


*****aaaB 


*  *\ 

•  ao  ” 


aaas9aaa««9eo3Q33 

o 


3i 
■3  * 

a 

c*** 


♦  a  a«  a ♦ # 

*  *  •  a  3  * 

*  o 
oa 


*  «  a  « 

oeaaaDQ  oeeeoeeQ 
*  «  *  a  a  a 


*  SS3* 

*  O  *  fr 

CO  *  * 


a  « 

oa 

#03 

*  *  o 

*  a  3  3  #  # 

♦  ♦  *  * ooaa  * 


e>*m 

*  *  *  *  c  3  oa 

oae^naao  oanaonao  * 

ft  It  •  tt  *  * 

*  *  *  * 


Figure  Ha  -  True  <0)  and  estimated  (*)  boundaries  Figure  lib  -  True  (0)  and  estimated  (*)  boundaries  for 

for  perfect  ellipse  using  curvature  algorithm,  perfect  ellipse  using  Kalman  filter  algorithm. 
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Figure  13  -  Est  iin.it cd  Cloud  Itoiind.ni  v 


